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The velocity fluctuations in a spherical shell arising from sinusoidal perturbations of a Keplerian shear flow with a free 
, amplitude parameter s are studied numerically by means of fully 3D nonlinear simulations. The investigations are performed 
^ ' at high Reynolds numbers, i.e. 3000 < Re < 5000. We find Taylor-Proudman columns of large eddies parallel to the rotation 

axis for sufficiently strong perturbations. An instability sets in at critical amplitudes with Scrit oc Re~^. The whole fiow 

turns out to be almost axisymmetric and nonturbulent exhibiting, however, a very rich radial and latitudinal structure. The 
' Reynolds stress {u'^u'^) is positive in the entire computational domain, from its Gaussian radial profile a positive viscosity- 

alpha of about 10 * is derived. The kinetic energy of the turbulent state is dominated by the azimuthal component (u^ ) 
■ whereas the other components are smaller by two orders of magnitude. Our simulations reveal, however, that these structures 
' disappear as soon as the perturbations are switched off. We did not find an "effective" perturbation whose amplitude is 
-««.^ , such that the disturbance is sustained for large times (cf. Dauchot & Daviaud 1995) which is due to the effective violation 
' of the Rayleigh stability criterion. The fluctuations rapidly smooth the original profile towards to pure Kepler flow which, 

therefore, proves to be stable in that sense. 

O , Key words: Turbulence theory, nonlinear hydrodynamics, angular momentum transport 

• • 1. Introduction 
> 

• I— I , . 

^> , The reason for the massive transport of angular momentum which enables the formation of the protoplanetary disk 
^ is the key question in accretion disk physics. Mostly this mechanism operates at very high Reynolds numbers so 

5^ that the influence of molecular viscosity is assumed to be negligible. This raises the problem of explaining efficient 
transport in nearly inviscid flows, which are still a challenge for numerical simulations because of the need for high 
numerical resolution due to the rich small-scale structure. 

It is widely accepted that small-scale turbulence emerging from shear flow instabilities can lead to enhanced 
outward transport - in opposition to the properties of convection phenomena in thin Keplerian disks. For the latter 
case several hydrodynamic computations have revealed an inwards transport of angular momentum, i.e. towards 
the high rotation rates. Ruden et al. (1988) estimated that a modest effective viscosity will result from the largest 
convective eddies. Ryu & Goodman (1992) started to find negative correlations {u'^u'^) for convection in Keplerian 
disks, and Cabot & Pollack (1992) found at low Reynolds numbers that the Reynolds stress can change sign with 
increasing rotation rate. Kley et al. (1993) derived from their numerical studies that large-scale convective motions 
can lead to an inward flux of angular momentum. Negative values for the mentioned cross correlation also appear 
in the simulations of Stone & Balbus (1996) probing the role of vertical convective motions to provide angular 
momentum transport in a Keplerian disk. In Riidiger et al. (2001) it is shown that an anisotropy in the form of 
dominating radial velocity fluctuations 

{<) > «) (1) 

should be responsible for negative cross correlations and v. v. Anisotropic turbulence fields under the influence of 
global rotation exhibits extra terms in the cross correlations which do not vanish for rigid rotation ('A- 

effect', see Riidiger 1989). Its sign depends on anisotropics in the turbulence field such as in (|l|). The question 
arises whether pure Keplerian shear fiow instabilities also fit this concept, although they are known to be linearly 
stable according to the Rayleigh criterion. 

For the solar nebula differential rotation alone has been suggested by DubruUe (1993) as a possible explanation 
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perturbations to the mean shear flow, which can occur due to pressure fluctuations or material impact, could lead 
to an instability independent of disk properties like temperature, pressure or density. 

Such finite amplitude perturbations for shear flows have been investigated numerically mainly in plane Couette 
flow, which is known to be linearly stable and due to the lack of curvature of the velocity profile the equations 
which have to be simulated are simplified significantly. Orszag & Kells (1980) have shown that a transition to 
turbulence at Re — 1250 requires fully 3D disturbances, but since they have used a rather low numerical resolution 
(16 X 33 X 16) they could not simulate the fully developed turbulent state. Finite-amplitude steady-state solutions 
have been found by Nagata (1990) for Re = 125. Lundbladh & Johansson (1991) studied the development of 
turbulent spots in plane Couette flow for Reynolds numbers between 300 and 1500 by means of direct numerical 
simulations, and actually they found turbulence to be sustained for sufficiently high Reynolds numbers. Re > 375. 
Corresponding experiments have been performed by Dauchot & Daviaud (1995) for Re < 400 to determine the 
critical amplitude e and found a power-law behavior 

e - (Re - Rcnl)"" with 0.3 < a < 0.8, (2) 

where Rcnl is the nonlinear critical Reynolds number, below which no sustained spots have been observed. This 
expression, valid only in the neighborhood of Rcnl, resembles a transition in terms of critical phenomena. However, 
for large Reynolds numbers an asymptotic behavior like 

e - Re"^ (3) 

holds. The parameter (3 was shown by Dubrulle & Zahn (1991) to be 1/3 for plane Couette flow. 

As has been pointed out by Balbus & Hawley (1996) the interaction between the velocity fluctuations and the 
background mean flow differs fundamentally for cartesian shear-flows and differential rotation. They used inviscid 
simulations of the Euler equations at moderate resolution of 64'^, but while turbulence was sustained in shear layers 
they found no numerical evidence that Keplerian flow at high Reynolds numbers is nonlinearly unstable. 

For accretion disks this problem now seems to have been solved by the finding that even a pure Keplerian shear 
fiow can be unstable when weak magnetic fields are present (Balbus & Hawley 1991). This magnetic shear flow 
instability does only operate in magnetically well coupled disks if the magnetic field is neither too weak nor too 
strong. 

In this paper we present results concerning the nonlinear stability of a perturbed Kepler flow for Reynolds 
numbers in the range 3000 < Re < 5000 obtained by fully nonlinear simulations of the Navier-Stokes equation. 
First we started with some sinusoidal internal noise perturbing the mean Keplerian shear flow and determined the 
critical amplitude for a given Reynolds number. In some sense, the modulation of the Kepler flow is used in order 
to simulate force-driven turbulence. Afterwards, when the kinetic energy has equilibrated, the perturbations have 
been switched off so that a pure Kepler flow remains. The kinetic energy decays very fast within a small fraction of 
a viscous time, thus our simulations finally confirm that a Kepler flow is nonlinearly stable even at high Reynolds 
numbers. 



2. Mathematical formulation 

We focus on the local and global properties of turbulence arising from finite amplitude perturbations. We consider 
a differentially rotating spherical shell of an incompressible fluid with inner radius Ri and outer radius i?o, which 
enables us to study the large scale structure as well. 

The mean flow is assumed to be Keplerian for large radii but perturbed by sinusoidal variations, i.e. 

n = no^^±^^^tt2^e^, (4) 

thus we are considering the noise induced by the perturbations as being part of the equilibrium profile (see Fig. |l]). 

This rotation law depends only on the distance zu = rsinO to the rotation axis, where the parameter wq is given 
by tuq = (i?o — Ri)/'2.. The perturbations we use are characterized by the parameter e describing the amplitude of 
the fluctuations and the number of periods 5 /it placed in a half shell diameter. This choice of localized perturbations 
allows for a smooth transition to the unperturbed Keplerian shear flow by decreasing the amplitude e. Because 
of numerical requirements we have restricted ourselves to (5 = 127r which can give rise to small-scale structures of 
radial dimensions of 1/24. 

According to the Rayleigh criterion a purely Keplerian shear flow is linearly stable with respect to infinitesimal 
axisymmetric perturbations, however due to fluctuations this stability criterion can be violated locally. It happens 
locally for all e values exceeding 0.007.|^ 

^As for finite value of e the basic flow may have a number of inflection points it is also worth to mention Rayleigh's inflection point 
theorem after which it is a necessary condition for instability to infinitesimal disturbances (for parallel shear flows) that an inflection 
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Fig. 1; LEFT: The mean Keplerian shear flow (dotted lino, e = 0) is perturbed by sine- like fluctuations. The full line 
in the diagram represents the resulting velocity proflle for the values e = 0.125 and 5 = 12it, the latter corresponding to 
approximately 12 periods. RIGHT: Negative values indicate local violations of the Rayleigh stability criterion 



Since this instability docs not depend on the behavior of temperature, pressure and density we do not have 
to take into account the temperature equation. Moreover, for simpUcity we assume the density to be constant 
throughout the whole shell. 

Length and time are normalized with respect to the difference in radii and to the viscous diffusion time r = 
(i?o - Rif/iy, i.e. 

r = {R,-R,)r, t = i^^^Mf (5) 

V 

Then the velocity field is normalized by the rate of the prescribed differential rotation, 

U = (iio - iii)OoM. (6) 

Expressed in terms of these variables the Navier-Stokes equation for the fluctuations u' =u — uq becomes 
du' 

— — Aw' = - gradp + Re {u' x rot u' + uo x rot u' + u' x rot uq) (7) 

(the hats are now omitted) with the Reynolds number 

^^^ Oo(i?o-i?0^ (8) 

Since the vector u is divergence- free it can be represented by toroidal and poloidal components (Chandrasekhar 
1961) 

u' = r X grad^^^ + rot^r x grad(— ) (9) 

We consider stress- free boundary conditions for the flow. Thus the cross-components Wrtj,, it re of the viscous 
stress-tensor -Kij = —py{u^ j + Uj ^) have to vanish, which can be expressed by the scalar potentials ^ and $ as 



r dr d'^r 

No normal flow is allowed at r = n, hence the potential ^' vanishes at the boundary, 
fl/. — n 
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3. Numerics 



The numerical simulations have been performed by the hydrodynamical part of the spectral magneto-convection 
code which is described in full detail in Hollerbach (2000). The spectral expansion uses spherical harmonics on the 
2-sphere and Chebyshev polynomials in the radial direction 

^ = ^ ^ik,l,m)Tk-iix)Pl"'^icos9)e'""*' (12) 

k,l,rn 

with the radius r mapped to the interval [— 1 , 1] via 

r = ^- + ^-x. (13) 

The radius ratio has been fixed to 1/5 throughout while setting the shell diameter to unity: Tq — rj = 1. In all 
diagrams with a radial dependence we show plots versus the ratio r/vo, normalizing the outer radius to unity. 

The diffusive terms in combination with the boundary conditions are treated implicitly in spectral space, whereas 
the nonlinear terms are evaluated explicitly on a grid of collocation points. The timestepping is performed by a 
modified second order Runge-Kutta scheme: first a predictor step calculates an estimated value, and afterwards 
the resulting spectral coefficients are used for evaluating again the velocity field, which yields a corrected value in 
a second step. 

Because of memory requirements we had to run these simulations on a Convex HPP-1200, where a single 
time step for a numerical resolution of 40 Chebyshev-, 70 Legendre polynomials and 3 Fourier modes required 
approximately 5 seconds of CPU time. At Reynolds numbers of Re = 3000 — 5000 the code runs stably with a 
time-step of 3 x 10~^. Thus a run simulating the flow over a period of half a viscous time can be performed within 
one day of CPU time. In order to avoid aliasing effects when transforming from real space to spectral space and 
vv. we used a collocation grid of the following dimension: dim(ri, 9j, 4>k) = 80 x 105 x 5. 



4. Results and discussion 



We have run simulations in the range Re= 3000 — 5000 with different values for the amplitude e so that the 
corresponding minimum perturbation can be determined. A single run where the instability leads to a nondecaying 
state is discussed in full detail in the following subsection. 



4.1. The exemplary case He— 3500 

Figure |2| shows the temporal evolution of the total kinetic energy for a run with Re= 3500 and several amplitudes 
e = 0.1, 0.125, 0.15 starting from arbitrary initial fields. For e — 0.15 the solution grows by one order of 
magnitude during a time of about 10% of a viscous time after which a turbulent state has been reached. The 
largest contribution to the kinetic energy is due to the azimuthal component, whereas the radial and meridional 
components yield contributions which are both smaller by almost two orders of magnitude. A perturbation with 
e — 0.125 still leads to a turbulent configuration, but it is close to the critical amplitude since the kinetic energy 
decays for e = 0.1. 

Figure |^ demonstrates the dominant role of the azimuthal velocity fiuctuations in the small-scale flow-pattern. 
The turbulence field is highly anisotropic with a massive "Austausch" in the ^-direction. As already stressed by 
Biermann (1951) such a flow field contributes to a positive correlation (m^m^) under the influence of basic rotation 
(cf. Riidiger 1989, p. 254) - in contrast to convective patterns, where negative correlations do appear. 

As the next step we discuss the spectral distribution of energy. The resulting flow pattern is nearly axisymmetric 
and steady, higher modes are suppressed by several orders of magnitude: Eq : Ei : E2 = I ■ 10"^ : 10^"*. 
Simulations with higher azimuthal resolution lead to the same quantitative result. For stronger amplitudes e the 
non-axisymmetric modes become much stronger, but in addition to the numerical difficulties that we encountered 
in simulating those flows such strong perturbations seem to be physically irrelevant. 

The Chebyshev spectrum (Fig. |4|) behaves like for small k and shows strong peaks at fci = 1 and k2 — 22, 
the former corresponding to the global torus-like structure of the flow. The latter peak determines the radial 
dimension of the turbulence elements Icoir — 1/^2 = 0.045. 

As can be seen from the energy spectrum for the Legendre polynomials the fiow is symmetric with respect 
to the equatorial plane. In general, even Z-modes are stronger by one order of magnitude with a maximum at 
Z = 46 in comparison with odd modes resulting in a total parity of 0.98. The flow is mainly conflned to the region 
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Fig. 2: The kinetic energy reaches a turbulent equilibrium for e = 0.15, 0.125, but decays for e = 0.1. The case e = 0.125 is 
close to the critical amplitude. Time is measured in units of the viscous time. 
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Fig. 3: The kinetic energy of the turbulent state is dominated by the azimuthal component, whereas the other components 
are smaller by two orders. 
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Fig. 4: The Chebyshev spectrum of the kinetic energy for Re= 3500 and e = 0.15 decreases like for large modes. The 
peak at fc = 22 corresponds to the radial dimension £corr = 0.045 of a single eddy. 
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Fig. 5: Left: Contour plot of the azimuthal velocity component. Right: Meridional structure of the velocity field. Re= 3500, 
e = 0.15. The eddies form Taylor-Proudman columns in the range 0.4 < m < 0.8. 
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0.0008 I 1 1 r-— 




RADIUS 

Fig. 6: The cross correlation of the velocity field has a Gauss-type structure in radius. The dotted lines show analytical 
curves exp— (ro — t47o)^/(2(7^) used for fitting the minima and maxima: vjo ~ 0.68, a ~ 0.204. e — 0.15, Re = 3500. 

numbers can also be seen from the contour plot of the radial flow component. We find 9 eddies in the radial and 
up to 18 eddies in the vertical direction depending on the height of the shell (Fig. ||). 

Despite the fact that the total flow looks rather inhomogeneous from a global point of view we make some 
rough estimates for its eddy viscosity which will agree approximately with the viscosity derived later on. From 
dimensional analysis one is led to 

= UT • 4orr, (14) 

V 

where €corr and ut denote a characteristic radial dimension and velocity. If we insert in this relation a typical eddy 
velocity and the shell diameter we end up with a rather small value oi v^/v k, 0.03. 

Now we turn to the angular momentum transport which is regulated by the dissipative stress. Figure ^ shows 
the Reynolds stress {u'ru' ^) as a function of the radius calculated by integrating over the 2-sphere for the amplitude 
£ = 0.15. Strong correlation of the velocity field can be found only in the region 0.5 < r < 0.9, the fluctuations 
being exponentially damped outside, but they are positive in the entire computational domain. Indeed after Fig. 

there is no exception from the dominance of {u'^) in striking difference to relation (|l|) valid for convective Kepler 
disks. 

The maxima and minima arc forming a Gaussian, 
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Fig. 7; The intensity ratio in the small-scale flow pattern of the model corresponding to Fig. 6. The azimuthal perturbations 
dominate the radial ones everywhere in opposition to the condition (1) valid for convection. 
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Fig. 8: The viscosity agg is positive and nearly constant between 0.5 < r < 0.9. The value of agg differs by one order for 
the amplitudes e = 0.15 and e = 0.125. 

with a standard deviation of about a ~ 0.204. We proceed to derive the eddy viscosity i/t by integrating ( p^ with 
respect to the radius and normalizing it with the unperturbed Keplerian shear flow Oxep: 



v \ dm 



cp 



-Re 



(16) 



This prescription yields in the case Re — 3500 for the ratio of turbulent and molecular viscosity vt/i^ — 0.15. This 
turbulent viscosity will now be related to the characteristic length and velocity. Shakura & Sunyaev (1973) have 
introduced a viscosity alpha using the scale height H{m) = -y/r^ — m'^. All uncertainties about the turbulence has 
been put into this single parameter which allows to construct a variety of a-type models in accretion disk theory, 

j/T -aggf^Kcpi?^- (17) 

The resulting radial structure for agg is shown in Fig. [s] for e — 0.15 and e — 0.125. It is positive everywhere. 
A complete volume integration yields the value agg = 3.75 x 10^'*. In the "turbulent" region 0.5 < r < 0.9 the 
viscosity alpha is nearly constant, while the value for e = 0.125 is by one order smaller than for e = 0.15. Those 
low values for all viscosity parameters indicate that we have not found a truly turbulent state, which is confirmed 
by the fact that the velocity field is stationary. In the next section we study the question whether this state can 
be sustained for a purely Keplerian rotation law. 



4.2. Perturbations as initial fields 



After the Rayleigh criterion a pure Keplerian shear fiow is stable against axisymmetric, infinitesimal perturbations. 
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Fig. 9: The kinetic energy {u' ) and the enstrophy (tu' ) decay exponentially within t — 0.1 r as soon as the perturbations 
£ are switched off. 



the nonlinear turbulent state discussed in the foregoing subsection can give rise to such a nonlinear instability the 
resulting flow was used as initial field for a simulation when the perturbations have been switched off, i.e. s = 0. 
Notice that the initial state contains weak non-axisymmctric contributions. 

It turns out that the kinetic energy (m'^) and the enstrophy (o;'^) decay exponentially very fast within a time of 
0.15 diffusion times, r, corresponding to approximately 525 orbits (Fig. ||). Within a few orbits the eddy viscosity 
decreases strongly since the rotational support from the sinusoidal perturbations is missing. After this short period 
it has reached approximately the value predicted by our estimate of the foregoing section using dimensional analysis: 
i/t/}^ ~ 0.02. Then it decays like i^t/v ~ 6"^^*/^. 

Both the enstrophy and the energy dissipation rate d{u''^)/dt are used for a comparison with the decay of 
homogeneous turbulence by rewriting the energy balance equation 

^^-i^. (18) 

where i^^q accounts for the effective nonlinear contributions stemming from the energy flux and for the energy 
injection by the forced differential rotation. As long as the flow shows strong local inhomogeneities this value 
remains nearly constant, but after t — 0.1 r, when the flow consists of several almost homogeneous domains, an 
exponential decay like i^^g/i' e~'^^'/'^ sets in. Hence, the turbulence does not survive the transition to a pure, 
undisturbed Kepler flow. 



5. Summary and conclusions 



The results for runs with several Reynolds numbers and with different perturbation amplitudes are summarized in 
Table 1 starting with Re= 3000. 

Since for stress-free boundary conditions the angular momentum is conserved, the kinetic energy does not 
decay to zero completely even for a purely Keplerian shear flow if the initial fields contain nonvanishing angular 
momentum. Therefore we started our series of runs at Re = 3000 without any perturbation. A calculation of the 
angular momentum yielded Lz ~ —5.4 x 10~^ in this case. We are considering the corresponding kinetic energy 
^kin ^ —5.89 x 10~^ as an offset for the subsequent runs with perturbations to the velocity profile e ^ rather 
than performing a transformation such that Lz = 0. 

The minimum value of the perturbation amplitude for reaching a turbulent state at Re= 3000 turned out to be 
£ = 0.15. The corresponding kinetic energy exceeded the energy offset due to initial angular momentum just a little 
bit, thus indicating that this case is close to neutral stability. The eddy viscosity accepted already a steady finite 
value of order 10^^ while the viscosity alpha reached approximately 10~^. Increasing the amplitude to e = 0.175 
all values grow by approximately one order of magnitude. Even for strong perturbations we never find a viscosity 
alpha larger than 10~^ indicating, therefore, no efficient mechanism of angular momentum transport in this model. 
The kinetic energy has been vastly dominated in all runs by the azimuthal component throughout. 

The onset of turbulence is mainly governed by the Reynolds number, for which one usually assumes that the 
critical amplitude behaves according to a power law like Ecrit ^ Re"''. For plane Couette flow an analytical 
treatment by DubruUe & Zahn (1991) yielded for this parameter /3 — 1/3. 

Using the critical amplitudes which we have determined in our runs we flnd roughly the following scaling 
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Table 1: The total kinetic energy, the eddy viscosity and the viscosity alpha for several combinations of Reynolds numbers 
and perturbation amplitudes. 



Via 


£ 




vt/v 


ass 


or\r\r\ 
oUUU 


A A 
U.U 


5.89 


xlO 


-4 


« 10-'"^ 




» 10-8 






0.125 


5.97 


xlO 


-4 


w 10-'"^ 




« 10-8 






0.15 


7.34 


xlO 


-4 


9.45 xlO" 


-3 


1.25 xlO- 


-5 




0.175 


8.55 


XlO 


-3 


1.03 xlO" 


-1 


1.84 xlO- 


-4 


3500 


0.1 


5.95 


xlO 


-4 


« 10-'^ 




« 10-8 






0.125 


9.37 


xlO- 


-4 


1.72 xlO" 


-2 


2.75 xlO- 


-5 




0.15 


1.56 


XlO 


-2 


1.50 xlO" 


-1 


3.75 xlO- 


-4 
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0.1 


5.85 


xlO 


-4 


w 10-'"^ 
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0.11 


1.73 


XlO 


-3 


3.70 xlO" 


-2 


6.63 xlO" 
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5.74 


XlO 


-4 
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0.085 


1.72 


xlO 


-3 


4.2 xlO" 


2 


7.52 xlO" 


-5 




0.09 


3.84 


xlO- 


-3 


7.71 xlO- 


-2 


1.53 xlO" 


'4 



This behaviour predicts that the instability will set in at high Reynolds numbers already for weak perturbations, but 
as our simulations have shown the nonaxisymmetric modes (being necessary for the occurrence of true turbulence) 
will be excited only for fairly strong perturbations. Accordingly, the turbulence was not sustained when the 
perturbations have been switched off, instead the eddy structures decayed very fast. Thus, our studies give 
no indications that Kepler rotation is nonlinearly unstable against finite amplitude perturbations (for Reynolds 
numbers up to 10^). 
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